TK!2
=v
#1
:n
L
:v
1000
:s
i
:u
ft
:d
ft
:c
length of pipeline
#2
:n
Q
:v
20
:s
i
:u
gpm
:d
gpm
:c
fluid rate
#4
:n
X1
:v
.3
:s
i
:c
/ input for bracketing
#5
:n
X2
:v
.4
:s
i
:c
\ procedure
#7
:n
x1
:v
.8236068025156
:s
o
:c
 /
#8
:n
x2
:v
1.24721360910153
:s
o
:c
|   triplet of values of
#9
:n
x3
:v
1.93262382478898
:s
o
:c
|   unknown  x  and objective
#10
:n
f1
:v
1109.17183465195
:s
o
:c
|   function  f  returned from
#11
:n
f2
:v
1013.82192572988
:s
o
:c
|   bracketing procedure
#12
:n
f3
:v
1248.89590474484
:s
o
:c
 \
#14
:n
xmin
:v
1.11726774189715
:s
o
:c
diameter giving minimum cost
#15
:n
costmin
:v
1003.00058665544
:s
o
:c
minimum annual cost
=u
=l
#1
:n
D_2nd
:d
in
:u
in
:v
.25
.287044475701847
.329578124123792
.378414319367648
.43448695960383
.498868326075083
.572789588409921
.657664348370408
.755115672263122
.867007128956061
.9954784270438231
1.14298633265317
1.31235165036323
1.50681316565986
1.7300895804699
1.98645062617284
2.28079871398955
2.6187626841543
3.00680544664253
3.45234757187584
3.96390919483858
4.55127294624867
5.22567102652779
6
E
#2
:n
cost_2nd
:v
21205.4373495424
12595.8937891136
7728.88473512323
4934.24895009567
3302.63727668545
2334.56548145275
1752.77070836336
1401.48029361543
1192.06099357184
1073.50742622561
1016.29662149917
1003.48090026427
1025.73468169597
1078.58724777611
1160.88510705895
1273.96413909939
1421.24817681271
1608.12016109254
1841.98410303459
2132.47737898075
2491.81767143739
2935.28494083424
3481.85048114449
4154.97484179372
E
#4
:n
k
:v
17
E
#5
:n
D
:v
.3
.4
.5618034
.8236068025156
1.24721360910153
1.93262382478898
1.93262382478898
1.50901702885931
1.24721362634371
1.08541022227338
1.09656424637357
1.12127966450115
1.11696724770761
1.11732003728717
1.11726885926601
1.11726774189715
1.11726662452941
E
#6
:n
cost
:v
10734.6219999275
4174.12413207817
1817.75074299717
1109.17183465195
1013.82192572988
1248.89590474484
1248.89590474484
1079.30328402401
1013.82192839229
1003.7965082141
1003.3316110659
1003.01260736534
1003.00065469391
1003.0005886743
1003.00058665564
1003.00058665544
1003.00058665712
E
=f
#1
:n
bracket
:c
Bracketing a minimum
:t
Procedure
:i
FUN,X1,X2
:o
x1,x2,x3,f1,f2,f3
#2
:n
Brent
:c
Inverse parabolic interpolation (Brent's method)
:t
Procedure
:i
FUN,a,b,c,tol
:o
x,fx
#3
:n
cost
:c
Objective Function
:t
Rule
:p
L,Q
:a
D
:v
cost
#4
:n
trace
:c
Optimization trace
:t
Procedure
:i
D,cost,Flag
#1
:n
bracket
:z
" Notation:  FUN    name of the objective function
"            X1,X2  initial search interval

" Description: The procedure searches for (and returns) set of three values
" of abscissa  x ,  x1 < x2 < x3  or  x1 > x2 > x3 , such that the corres- 
" ponding function values  f1 > f2  and  f3 > f2 .  This means that the
" minimum should be somewhere between  x1  and  x3 ,  and that an optimiza-
" tion procedure, such as golden section search 'Golden' or inverse parabolic
" interpolation 'Brent' can converge to it.

x1:= X1
x2:= X2
gold:= 1.618034                    " interval magnification ratio
glimit:= 100                       " maximum allowed magnification
f1:= apply(FUN,x1)
f2:= apply(FUN,x2)
if f1 > f2 goto A
  (x1,x2):= (x2,x1)                " securing downhill move from x1 to x2
  (f1,f2):= (f2,f1)
A: x3:= x2 + gold*(x2-x1)
   f3:= apply(FUN,x3)
B: if f2 < f3 return
   xf1:= (x2-x1)*(f2-f3)           " finding x4 by parabolic extrapolation
   xf3:= (x2-x3)*(f2-f1)
   if f1=f3 then f1:= .999*f1               " preventing division by 0
   x4:= x2 - ((x2-x3)*f3 - (x2-x1)*f1) / (f3-f1)
   xlimit:= x2 + glimit*(x3-x2)             " current step limit
   if (x2-x4)*(x4-x3)<=0 goto E
     f4:= apply(FUN,X4)
     if x4 >= x3 goto C
       (x1,f1):= (x2,f2)                    " minimum between x2 and x3
       (x2,f2):= (x4,f4)
       goto B
     C: if f4 <= f2 goto D
          (x3,f3):= (x4,f4)                 " minimum between x1 and x4
          goto B
        D: x4:= x3 + gold*(x3-x2)
           f4:= apply(FUN,x4)
           goto H
   E: if (x3-x4)*(x4-xlimit) <= 0 goto F
        f4:= apply(FUN,x4)                  " x4 is between x3 and xlimit
        if f4 >= f3 goto H
          (x2,f2):= (x3,f3)
          (x3,f3):= (x4,f4)
          x4:= x3 + gold*(x3-x2)
          f4:= apply(FUN,x4)
          goto H
      F: if (x4-xlimit)*(xlimit-x3) < 0 goto G
           x4:= xlimit
           f4:= apply(FUN,x4)
           goto H
         G: x4:= x3 + gold*(x3-x2)
            f4:= apply(FUN,x4)
H: (x1,f1):= (x2,f2)                        " eliminating x1
   (x2,f2):= (x3,f3)
   (x3,f3):= (x4,f4)
   goto B
/E
E
#2
:n
Brent
:z
" Notation:  FUN    name of the objective function
"            a,b,c  values of  x  bracketing the minimum
"            tol    desired relative precision of  x  giving minimum
"                   FUN (shouldn't be less than 1e-8)

" Description: The procedure returns argument  x  and function value  fx 
" corresponding to the minimum of the objective function FUN within the
" interval of the arguments values between  a  and  c .  It uses Brent's 
" algorithm of repeated inverse parabolic interpolation, i.e fitting a para-
" bola through three bracketing points and finding its minimum; when the
" parabolic interpolation fails to indicate convergence to minimum, the
" algorithm switches to the golden section search.

" Note: Maximum number of iteration steps is set in the statement labeled E.

gold:= .381966
eps:= 1e-10
(a,b):= (min(a,c),max(a,c))
(v,fv):= (b,apply(FUN,b))
(w,fw):= (v,fv)
(x,fx):= (v,fv)
e:= 0
iter:= 0
LOOP:
iter:= iter + 1
m:= (a+b)/2
tol1:= tol*abs(x) + eps
tol2:= 2*tol1
if abs(x-m) <= (tol2-(b-a)/2) return
if abs(e) <= tol1 goto A

  r:= (x-w)*(fx-fv)                      " trying parabolic fit
  q:= (x-v)*(fx-fw)
  p:= (x-v)*q - (x-w)*r
  q:= 2*(q-r)
  if q > 0 then p:= -p
  q:= abs(q)
  etemp:= e
  e:= d
  if or(abs(p)>=abs(.5*q*etemp),p<=q*(a-x),p>=q*(b-x)) goto A

    d:= p/q                               " taking parabolic step
    u:= x + d
    if or((u-a)<tol2,(b-u)<tol2) then d:= abs(tol1)*sgn(m-x)
    goto B
A: if x < m then  e:= b-x  else  e:= a-x

   d:= gold*e                             " taking golden section step
B: if abs(d) < tol1 then  u:= x + abs(tol1)*sgn(d)  else  u:= x + d

   fu:= apply(FUN,u)                      " function evaluation at new point

   if fu > fx goto C                      " incorporating new point
     if u < x then  b:= x  else  a:= x    " into the search scheme
     (v,fv):= (w,fw)
     (w,fw):= (x,fx)
     (x,fx):= (u,fu)
     goto E
   C: if u < x then  a:= u  else  b:= u
      if and(fu>fw,w<>x) goto D
        (v,fv):= (w,fw)
        (w,fw):= (u,fu)
        goto E
      D: if or(fu<=fv,v=x,v=w) then (v,fv):= (u,fu)
E: if iter<50 goto LOOP
error:= convergence_failure
/E
E
#3
:n
cost
:z
cost = .45*L + .245*L*D^1.5 + 325*sqrt(hp) + 61.6*hp^.925 + 102
hp = L * (4.4e-8*Q^3/D^5 + 1.92e-9*Q^2.68/D^4.68)
call trace(D,cost,'record)
/E
E
#4
:n
trace
:z
" Notation:  D     unknown variable value
"            cost  objective function value
"            flag  initialization/recording flag

" Description: This is an example of a procedure for tracing the optimiza-
" tion process.  Here are the steps for incorporating it into your model:

" 1. Select the variable you want to record and make them Input Variables
"    in the  trace  function.

" 2. Add a variable  Flag  as an input variable.

" 3. Select the list names under which you want to record the trace infor-
"    mation.  Blank these lists in the initialization part of this procedure
"    and use them on the left hand side of the assignment statement in the
"    recording part. (We use the same names for the variables and lists
"    in this example.)

" 4. Set up some list element as a record counter.  Set it to zero in the
"    initialization part and increment it in the recording part. (This 
"    example uses k[1].)

" 5. Include an initialization call of the tha function  trace  before the
"    optimization procedure call on the Rule Sheet of Function Subsheet.
"    All the arguments may be set to zero, e.g.  call trace(0,0,0) .

" 6. Include the recording call(s) in the Function Subsheet evaluating the
"    objective function.  The calls must have as a last argument a symbolic
"    constant  'record .

" 7. Set up an Interactive Table for displaying the trace lists.

" Note: You may want to display more information than just the values of 
" unknowns and the objective function.  Any such information can be sent
" to the  trace  via input variable arguments. 

if Flag='record goto RECORD
'k[1]:= 0                         " initialization
  call blank('cost)
  call blank('D)
  return
RECORD:
'k[1]:= 'k[1] + 1                 " recording the objective function call
'cost['k[1]]:= cost
'D['k[1]]:= D
/E
E
=r
#1
:r
" ***  Single-variable Optimization, Part 4  ***
:s
C
#3
:r
" This model demonstrates the use of bracketing the minimum and finding it
:s
C
#4
:r
" by Brent's inverse parabolic interpolation method:
:s
C
#6
:r
call trace(0,0,0)
#7
:r
call bracket('cost,X1,X2;x1,x2,x3,f1,f2,f3)
#8
:r
(xmin,costmin) = Brent('cost,x1,x2,x3,1e-6)
#10
:r
" Notes:
:s
C
#12
:r
" 1. The problem deals with the design of a pipeline for a chemical plant
:s
C
#13
:r
"    satisfying the condition of minimum annual cost of pipe, pump and 
:s
C
#14
:r
"    pumping (see the description in model OPTIM1A and alternative methods
:s
C
#15
:r
"    of solution in models OPTIM1A and OPTIM1B).  The objective function
:s
C
#16
:r
"    is defined in the rule function  cost .
:s
C
#18
:r
" 2. The function  bracket  takes two guess points (in our example  X1=.3 
:s
C
#19
:r
"    and X2=.4, see the Variable Sheet), searches the objective function, 
:s
C
#20
:r
"    and generates three points with a minimum value of  f  somewhere between
:s
C
#21
:r
"    the smallest and the largest  x  (check the description in the Function
:s
C
#22
:r
"    Subsheet, and the values  x  and  f  in the Variable Sheet).  Press F7
:s
C
#23
:r
"    to display the graph of the objective function from model OPTIM1A; you
:s
C
#24
:r
"    can see that we have chosen the values of  X1  and  X2  deliberately so
:s
C
#25
:r
"    that the function  bracket  had to "work" to locate an interval with the
:s
C
#26
:r
"    minimum point inside, and it did!
:s
C
#28
:r
" 3. Next, the output from the function  bracket  is fed into the function
:s
C
#29
:r
"    Brent which finds, with desired precision of 1e-6, a value of  x
:s
C
#30
:r
"    for minimum  cost .  The function  Brent  is described in the Function
:s
C
#31
:r
"    Subsheet.
:s
C
#33
:r
" 4. Also included in the model is function  trace  that records the calls to
:s
C
#34
:r
"    the function  cost  (see the description in the Function Subsheet).
:s
C
#35
:r
"    Press  = T > >  to see the trace records.  The first six entries show
:s
C
#36
:r
"    the work of the function  bracket , the rest reflects the optimization
:s
C
#37
:r
"    process controlled by  Brent .
:s
C
=p
#1
:n
second
:t
Line chart
:d
EGA
:i
""
:s
Yes
:z
None
:g
No
:l
None
:h
"Pipe Diameter"
:v
"Cost"
:x
D_2nd
:u
#1
:y
cost_2nd
:f
1
:p
Lines
:c
"*"
E
=t
#1
:n
trace
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
""
:u
#1
:l
D
:w
20
:h
""
#2
:l
cost
:w
20
:h
""
E
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tr,f,0,1,0,1,0
